function polynomial = computePolynomial(M, N)

% compute the unique terms of a M order polynomial with N variables

polynomial = ones(N^M, M);
for i = 2:size(polynomial, 1)
    polynomial(i, :) = polynomial(i-1, :);
    polynomial(i, end) = polynomial(i, end) + 1;
    for j = 1:M-1
        if(mod(polynomial(i, end-j+1), N+1) == 0)
            polynomial(i, end-j+1) = 1;
            polynomial(i, end-j) = polynomial(i, end-j) + 1;
        end
    end
end

for i = 1:size(polynomial, 1)
    polynomial(i, :) = sort(polynomial(i, :));
end

polynomial = unique(polynomial, 'rows');
